Initialisation

## [1] "number of cores available = 16"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[1] + eta)/(1+exp((phi[2]-t)/phi[3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,90,5),
              omega2 = c(0.5,0.1,0.01),
              #Survival data,
              nu2 = 0.5,
              a = 90,
              b = 50,
              alpha = 7,
              beta = 10)

F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(60,120, length.out = 10) #value of times

# --- longitudinal data --- #
G = 40 ; ng = 4 #nombre de groupe et d'individu par groupe
dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs
dt_NLME %>% ggplot(aes(t, obs, col = gen, group = id)) + 
  geom_point() + geom_line() + theme(legend.position = 'null')

# --- Survival data --- #
dt_SF <- SF_obs(dt_NLME, param, m)
dt_SF %>% ggplot(aes(obs, fill = factor(gen))) + 
  geom_histogram(position = 'dodge', bins = 20) + theme(legend.position = 'null')

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

sigma2_a <- sigma2_b <- sigma2_alpha <- 0.1

Phi <- function(sigma2, rho2, mu, omega2, nu2, a, b, alpha, beta)
  c(- 1/(2*sigma2), -1/(2*rho2), -1/(2*omega2), G*mu/omega2,  #1, 2, 3, 4
    -1/(2*nu2), a/sigma2_a,         -1/(2*sigma2_a),          #5, 6, 7
                b/sigma2_b        , -1/(2*sigma2_b),          #8, 9
                alpha/sigma2_alpha, -1/(2*sigma2_alpha)  )    #10, 11

zeta <- function(beta, eta, phi, gamma, a, b, alpha)
{
  lbd <- function(t,g) t^{b-1} * exp(alpha*m(t, eta[g], phi[g,]))
  P1 <- 1:nrow(dt_SF) %>% sapply(function(i) integrate(function(t) lbd(t,dt_SF$gen[i]) , 0, dt_SF$obs[i])$value )
  
  P2 <- 1:nrow(dt_SF) %>% sapply(function(i) dt_SF$U[i]*exp(beta*dt_SF$U[i] + gamma[dt_SF$gen[i]] )  )
  
  sum(dt_SF$U) -  b*a^-b * sum(P2*P1)
}

S <- function(eta, phi, gamma, a, b, alpha)
{
  beta <- uniroot(function(beta)zeta(beta, eta, phi, gamma, a,b,alpha), lower = -100, upper = 100)$root
  
  s <- c(mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ), #1
         mean(eta^2),                                               #2
         apply(phi^2, 2, mean), apply(phi  , 2, mean),              #3 et 4
         mean(gamma^2),                                             #5
         a, a^2,                                                    #6 et 7
         b, b^2,                                                    #8 et 9
         alpha^2, alpha,                                            #10 et 11
         beta )                                                     #12

  attr(s, '1') <- 1
  attr(s, '2') <- 2
  attr(s, '3') <- 2 + 1:ncol(phi)
  attr(s, '4') <- 2 +   ncol(phi) + 1:ncol(phi)
  attr(s, '5') <- 2 + 2*ncol(phi) + 1
  
  for(i in 6:12) attr(s, as.character(i)) <- attr(s, as.character(i-1)) + 1 
  return(s)
}
#=============================================================================#
loglik.phi <- function(x, eta, Phi)
{
  Phi[1] * sum((Y - get_obs(m, dt_NLME, eta = eta, phi = x) )^2) +
    Phi[3] * apply(x^2, 2, sum) + Phi[4] * apply(x  , 2, sum)
}
loglik.eta <- function(x, phi, Phi) 
{ 
  Phi[1] * sum((Y - get_obs(m, dt_NLME, eta = x, phi = phi) )^2) + 
    Phi[2] * sum(x^2)
}

loglik.gamma <- function(x, Phi) Phi[5]  * sum(x^2) + ng*sum(x) #Correction fait à la va vite
loglik.a     <- function(x, b, Phi) Phi[6]  * x + Phi[7] * x^2  + G*ng*log(b*x^-b)
loglik.b     <- function(x, Phi) Phi[8]  * x + Phi[9] * x^2  + (x-1)*sum(log(t))
loglik.alpha <- function(x, Phi) Phi[10] * x + Phi[11] * x^2

SAEM

initialisation

M <- 1 #nombre de simulation
u <- function(k) ifelse(k<200, 1, 1/(k-199) )

# ---  Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x + 2*runif(1) )

# --- Initialisation des chaines MC : Z_0 --- #list(attr(dt_NLME,'eta')),#
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, para$rho2)  %>% matrix(ncol = 1) ),
          phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ),
          #list(attr(dt_NLME,'phi')),#
          gamma = 1:M %>% lapply(function(i) matrix(rnorm(G, 0, para$nu2), ncol = 1) ),
          
          a = list(matrix(para$a)), #list(matrix(param$a)), #
          b = list(matrix(para$b)), #list(matrix(param$b)), #
          alpha = list(matrix(para$alpha)) )#list(matrix(param$alpha)) )#

Boucle

sim <- function(niter, Phih, eta, phi, gamma, a, b, alpha)
{
  M <- length(phi)

  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_para(niter, eta[[i]],   sd = .02,              loglik.eta, phi[[i]], Phih, cores = ncores-1 ))
  #                  |                                         |                     |
  phi <- 1:M %>% lapply( function(i) #                         |                     |
    MH_High_Dim_para(niter, phi[[i]],   sd = c(.03, .02, .02), loglik.phi, eta[[i]], Phih, cores = ncores-1 ))
  #                  |                                         |                     |
  gamma <- 1:M %>% lapply( function(i)#                        |                     |
    MH_High_Dim_para(niter, gamma[[i]], sd = .03,              loglik.gamma,         Phih, cores = ncores-1 ))

  a <- 1:M %>% lapply( function(i) MH_High_Dim_para(niter, a[[i]],     sd = .02, loglik.a,b[[i]],  Phih, cores = 1 ))
  b <- 1:M %>% lapply( function(i) MH_High_Dim_para(niter, b[[i]],     sd = .02, loglik.b,         Phih, cores = 1 ))
  alpha <- 1:M %>% lapply( function(i) MH_High_Dim_para(niter, alpha[[i]], sd = .02, loglik.alpha, Phih, cores = 1 ))

  list(eta = eta, phi = phi, gamma = gamma, a = a, b = b, alpha = alpha)
}

maxi <- function(S)
{
  list(sigma2 = S[attr(S,'1')],
       rho2 =   S[attr(S,'2')],
       mu =     S[attr(S,'4')],
       omega2 = S[attr(S,'3')] - S[attr(S,'4')]^2,
       nu2 =    S[attr(S,'5')],
       a =      S[attr(S,'6')],  b =      S[attr(S,'8')],
       alpha =  S[attr(S,'11')], beta =   S[attr(S,'12')] )
}
Estimation oracle
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23 nu2 a b alpha beta
0.0511813 0.1045811 4.932414 90.07442 4.980973 0.4093858 0.0762974 0.0081025 0.5462794 90 50 7 9.967896
niter <- 300
MH.iter <- 10
res <- SAEM(niter, MH.iter, u, para, Phi, S, Z, sim, maxi, eps = 1e-1, verbatim = T)
saveRDS(res, 'saem.rds')
gg <- plot(res)
## [1] "SAEM execution time = 00min 02sec"
Résultat de l’algo SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 nu2 a b alpha beta
Valeur réelle 0.0500 0.1000 5.0000 90.0000 5.0000 0.5000 0.1000 0.0100 0.5000 90.0000 50.0000 7.0000 10.0000
Valeur estimée 0.0517 3.4131 3.6631 89.9177 5.0576 1.7349 0.3362 0.2266 0.8383 107.1835 68.2192 23.7318 20.0341
Rrmse 0.0342 33.1310 0.2674 0.0009 0.0115 2.4697 2.3615 21.6611 0.6766 0.1909 0.3644 2.3903 1.0034